function [P,dP,hP] = logricepdf(x, A, sigma, selGrHess)
% [P] = LOGRICEPDF(X, A, SIGMA)
% Returns the logarithm of the rice probability density of the values X, 
% with amplitude A and a  noise level of sigma. Default: sigma= 1
%
% [P,dP,d2P] = LOGRICEPDF(X, A, SIGMA, selGrHess)
% Additionally compute the first (dP) and second (d2P) derivative of the
% log pdf with respect to the X (selGrHess(1) = true), A (selGrHess(2) = true) 
% and/or, sigma (selGrHess(3) = true). default: selGrHess = [true true true];
%
% Created by Dirk Poot, University of Antwerp, 13-8-2007

% 22-6-2011 : added support for expansion of sigma (bsxfun)
%             added some commments explaining the code.

if nargin<3
    sigma = 1;
    if nargin==1 && isequal(x,'printtables')
        print_c_tables();return;
    end;
end;

% see http://en.wikipedia.org/wiki/Rice_distribution
% pdf(x | A, sigma) = x/sigma^2 * exp( -(x^2 + A^2)/(2 sigma^2) ) * I_0(x*A/sigma^2) 
% with I_0 = modified bessel function of first kind with order 0.

rs2 = 1./sigma.^2;
bx = bsxfun(@times, rs2, x.*A);

% Compute bessel_0(bx)
% use  1  as third argument to scale bessel output by exp(-abs(real(bx))). This avoids overflow for large A,x . 
% Correct the logarithm afterwards by adding abs(real(bx)) .
usemex = true; % use much faster custom mex function?

if nargout>1
    % Compute gradient wrt x, A, and/or sigma
    if nargout>2
        if usemex
            [logbesX, besRat, besRat2m] = rice_besseli_funs_c(bx);
        else
            besX = besseli(0,bx,1);
            besRat = besseli(1,bx,1)./besX;
            besRat2 = besseli(2,bx,1)./besX;
            besRat2m = (besRat2 -2 .* besRat.^2 +1).*bx;
            logbesX = log(besX);
        end;
    else
        if usemex
            [logbesX,besRat]=rice_besseli_funs_c(bx);
        else
            besX = besseli(0,bx,1);
            besRat = besseli(1,bx,1)./besX;
            logbesX = log(besX);
        end;
    end;
    
    if nargin<4
        selGrHess = true(1,3);
    elseif ~islogical(selGrHess) || numel(selGrHess)~=3
        error('incorrect selGrHess, should be 3 element logical vector.');
    end
    % compute gradient.
    if selGrHess(1)
        dPdx = 1./x + bsxfun(@times, rs2 , A .* besRat - x);
    else
        dPdx = [];
    end
    if selGrHess(2)
        dPdA = bsxfun(@times, rs2  , x .* besRat - A);
    else
        dPdA = [];
    end;
    if selGrHess(3)
        dpds = bsxfun(@times,  bsxfun(@times, rs2, (A.^2 + x.^2)) - 2 * bx .* besRat - 2 , 1./ sigma);
    else
        dpds = [];
    end;
    dP = [dPdx(:) dPdA(:) dpds(:)];
    
    if nargout>2
        % also compute hessian.
        if selGrHess(1)
%             dPdxx = -1./x.^2 + rs2 .*( .5* rs2 .* A.^2 .*( 1 - 2* besRat.^2 + besRat2 ) - 1 );
            dPdxx =  (bsxfun(@times, rs2 ,  .5* A.*x .* besRat2m  - x.^2 )-1 )./x.^2;
            if selGrHess(2)
                dPdxA = bsxfun(@times, rs2 , (besRat +.5*besRat2m));
            end;
            if selGrHess(3)
                dPdxs = bsxfun(@times, rs2./sigma , (2*x - A.*(2*besRat +besRat2m)));
            end;
        end;
        if selGrHess(2)
%             dPdAA = bsxfun(@times, rs2 , bsxfun(@times, .5* rs2 , x.^2 .*( 1 - 2* besRat.^2 + besRat2 )) - 1 );
            dPdAA = bsxfun(@times, rs2 , .5*besRat2m.*x./A - 1 );
            if selGrHess(3)
                dPdAs = bsxfun(@times, rs2./sigma , (2 .* A - (2 .* besRat + besRat2m) .* x));
            end;
        end;
%         if selGrHess(2) && selGrHess(3)
%             dPdAs = bsxfun(@times, rs2./ sigma , 2*A - 2 .* x .* besRat + bsxfun(@times, rs2 , 2 .* A .* x.^2 .*besRat.^2 - A .* x.^2 .*(1 + besRat2)) ) ;
%         end;
        if selGrHess(3)
%             dPdss = bsxfun(@times, rs2 , 2 + bsxfun(@times, rs2 , (-3*(x.^2+A.^2) + bsxfun(@times, rs2 , x.^2.*A.^2 .*( 2 - 4 *besRat.^2 + 2 .* besRat2)) + 2*A.*x .* (3.*besRat))));
            dPdss = bsxfun(@times, rs2 , 2 + bsxfun(@times, rs2 , (-3*(x.^2+A.^2) + 2*A.*x .*( 3 *besRat + besRat2m)) ));
        end;
        if isequal(selGrHess,[true false false])
            hP = dPdxx(:);
        elseif isequal(selGrHess,[false true false])
            hP = dPdAA(:);
        elseif isequal(selGrHess,[false false true])
            hP = dPdss(:);
        elseif isequal(selGrHess,[false true true])
            hP = [dPdAA(:) dPdAs(:) dPdss(:)];
        elseif isequal(selGrHess,[true true true])
            hP = [dPdxx(:) dPdxA(:) dPdAA(:) dPdxs(:) dPdAs(:) dPdss(:)];
        else
            error('hessian with this selGrHess pattern not yet supported');
        end;
    end;
else
    % nargout<=1
    if usemex
        logbesX = rice_besseli_funs_c(bx);
    else
        logbesX = log(besseli(0,bx,1));
    end;
end;
% -(x^2+A^2)/(2*sigma^2) + x*A/sigma^2
% = - ( x - A)^2/(2*sigma^2)
P = bsxfun(@plus, log(x) ,log(rs2)) + bsxfun(@times, -.5.*rs2, (x-abs(A)).^2) + logbesX;


function print_c_tables()
cutoffbinidx=29;
limits = [0,41792,125506,203745,280270,356795,434113,513004,594255,678896,768745,867886,977719,1082396,1181955,1280388,1379346,1479994,1583388, ...
  1690698,1803477,1924490,2060206,2208398,2351328,2489750,2627784,2767519,2910366,3057509,(1099511627776/339765),(549755813888/160183),( ...
  1099511627776/301045),(1099511627776/280779),(1099511627776/259495),(274877906944/60497),(549755813888/112917),(1099511627776/210287),( ...
  549755813888/97421),(1099511627776/177993),(17179869184/2561),(1099511627776/150931),(274877906944/34509),(274877906944/30889),(1099511627776/111523), ...
  (1099511627776/99615),(137438953472/10713),(1099511627776/73577),(1099511627776/57665),(1099511627776/41865),(274877906944/5659),10485760000]'*( ...
  1/1048576);
origins = [0,167298,329251,484015,637065,790908,947117,1107259,1273151,1447641,1636631,1845605,2060115,2264351,2462343,2659734,2859340, ...
  3063382,3274086,3494175,3727967,3984696,4268604,4559726,4841078,5117534,5395303,5677885,5967875,699404,660131,621411,581824,540274,501483,467822, ...
  436121,405129,372835,341897,314835,288967,261592,235079,211138,185319,159281,131242,99530,64501,22636]'*(1/2097152);
funcs={[0,(-1),(1/4),0,(-0.156249999955908758E-1),(-0.354009593034203857E-8),0.173655525467850941E-2,(-0.178328120895867687E-4);(-0.781835663457938204E-1),(-0.960144744386072799E0),0.249404439710345575E0,( ...
  -0.496828187227357229E-2),(-0.154599061644540754E-1),0.824640767586818742E-3,0.169581019018079974E-2,(-0.140874769436868943E-3);(-0.150846398147436310E0),(-0.921741319201242991E0), ...
  0.247704910945054459E0,(-0.967926221102290059E-2),(-0.149925271570229497E-1),0.158764360990112645E-2,0.158500906586203519E-2,(-0.266931374785202992E-3);(-0.217523673089228052E0),( ...
  -0.885363435469714822E0),0.245079177304576671E0,(-0.140059817782421972E-1),(-0.142813248088064478E-1),0.225503806902989882E-2,0.141973103879107830E-2,(-0.369661556286934091E-3);( ...
  -0.280837967551910310E0),(-0.849837372315784601E0),0.241565683529992712E0,(-0.180442214754695628E-1),(-0.133503294248866953E-1),0.283236252300445354E-2,0.120976420802177423E-2,( ...
  -0.448269061908029931E-3);(-0.341887928584762494E0),(-0.814707758963825771E0),0.237175241648150075E0,(-0.218001307609726479E-1),(-0.122201745137864434E-1),0.331205750166953709E-2, ...
  0.964975023360291889E-3,(-0.500582126151962147E-3);(-0.401265940421560148E0),(-0.779757700436698066E0),0.231911120128883323E0,(-0.252498809164337899E-1),(-0.109137008257382334E-1), ...
  0.368382618808224143E-2,0.696383433444933857E-3,(-0.524969827921895259E-3);(-0.459468836803714640E0),(-0.744799971815103550E0),0.225761648078782005E0,(-0.283630505636834912E-1),( ...
  -0.945447314790304685E-2),0.393847623972151870E-2,0.415610823892325028E-3,(-0.521012186262270577E-3);(-0.516986824254234605E0),(-0.709633286911306517E0),0.218695545753311426E0,( ...
  -0.311047269328553097E-1),(-0.786666277889674771E-2),0.406830302983904597E-2,0.134688808813454018E-3,(-0.489574252000247823E-3);(-0.574534960016481360E0),(-0.673904040268208287E0), ...
  0.210628230643125622E0,(-0.334404721245479837E-1),(-0.616983633348759616E-2),0.406673984581501645E-2,(-0.134810751480013843E-3),(-0.432591065190680323E-3);(-0.633579788600083909E0),( ...
  -0.636772933209760478E0),0.201316459355650008E0,(-0.353371881500809467E-1),(-0.436447995467310965E-2),0.392421710149743690E-2,(-0.382967306429269271E-3),(-0.352076922291638023E-3);( ...
  -0.695068422760418838E0),(-0.597719943766344119E0),0.190530940883209826E0,(-0.366958936257076231E-1),(-0.247777161188102665E-2),0.362834709246623469E-2,(-0.593874618900729040E-3),( ...
  -0.252141808635828117E-3);(-0.754253080569088416E0),(-0.559902973371987951E0),0.179152731117026940E0,(-0.373436584361439481E-1),(-0.723901582622721301E-3),0.321564361825707702E-2,( ...
  -0.737525979133022546E-3),(-0.149899377489014057E-3);(-0.807115920935891688E0),(-0.526072368875704248E0),0.168229821871914824E0,(-0.373347216067160589E-1),0.732801431442176538E-3, ...
  0.276054528432908441E-2,(-0.809089911895653730E-3),(-0.619481645628286764E-4);(-0.855314249536425474E0),(-0.495302026418923216E0),0.157716962398132841E0,(-0.368256868507696533E-1), ...
  0.192641581713683869E-2,0.229506898951687493E-2,(-0.825682144118887061E-3),0.924201515361137883E-5;(-0.900567034735653565E0),(-0.466583774109562134E0),0.147439057152115133E0,(-0.359107961469697144E-1) ...
  ,0.289745085063618388E-2,0.183396208132727140E-2,(-0.800639466182675247E-3),0.640751411497053282E-4;(-0.943671312166470810E0),(-0.439482660620430282E0),0.137357465766390929E0,( ...
  -0.346551441815756307E-1),0.366366265182827549E-2,0.139151271225806972E-2,(-0.744107426866654470E-3),0.102971227539677895E-3;(-0.985162007095988377E0),(-0.413724407708870163E0),0.127462091999119241E0, ...
  (-0.331109669095426088E-1),0.423846448873131692E-2,0.979346141917129350E-3,(-0.664966270697137440E-3),0.127007124593464236E-3;(-0.102547598894415595E1),(-0.389096857689448399E0), ...
  0.117747622406787052E0,(-0.313217534575567617E-1),0.463438456963846562E-2,0.606354478274941376E-3,(-0.571108437085187437E-3),0.137826316500485158E-3;(-0.106504917305202173E1),(-0.365395590858980163E0) ...
  ,0.108198545845706628E0,(-0.293221340265484634E-1),0.486380891536483938E-2,0.278781993243306966E-3,(-0.469393620894544980E-3),0.137448053259690734E-3;(-0.110447891793012369E1),( ...
  -0.342337607384970054E0),0.987574960130841225E-1,(-0.271308734627363224E-1),0.493828654579741086E-2,0.636737674052007101E-7,(-0.365356774458055343E-3),0.128021390419002979E-3;(-0.114495584924498753E1) ...
  ,(-0.319341812093838572E0),0.892364587214447514E-1,(-0.247251438139326140E-1),0.486419640631579615E-2,(-0.229544325438311309E-3),(-0.262483249012788825E-3),0.111417544483324097E-3;( ...
  -0.118661193968042271E1),(-0.296492130339444821E0),0.797227137091353023E-1,(-0.221449507858403609E-1),0.464592637391539929E-2,(-0.402440156264094594E-3),(-0.166983649668371434E-3), ...
  0.899975722829580576E-4;(-0.122629157021705390E1),(-0.275589575383338434E0),0.710259458221373896E-1,(-0.196505651686442941E-1),0.432629452919068446E-2,(-0.507917948661560864E-3),( ...
  -0.901906460175337077E-4),0.682898299339372549E-4;(-0.126203211537531346E1),(-0.257552131686808948E0),0.635716023903795034E-1,(-0.174239546176567789E-1),0.396663971446936433E-2,( ...
  -0.557017227428638580E-3),(-0.350432135321145129E-4),0.496004998314259936E-4;(-0.129491780597628643E1),(-0.241664408189177789E0),0.570815976254179127E-1,(-0.154302602344914410E-1), ...
  0.359404794784155444E-2,(-0.568491208105170721E-3),0.342189528087561751E-5,0.342854978608623787E-4;(-0.132595977655200506E1),(-0.227322981275895507E0),0.513154967620122532E-1,( ...
  -0.136253519379656690E-1),0.322101246620665302E-2,(-0.554653418918603486E-3),0.293171936303944792E-4,0.220879600729745622E-4;(-0.135569114104207271E1),(-0.214205452540822898E0), ...
  0.461451144670446283E-1,(-0.119883184171841191E-1),0.285700596265420463E-2,(-0.523748249506133570E-3),0.455046565013440717E-4,0.127069374595914070E-4;(-0.138445939094395401E1),( ...
  -0.202102150042184346E0),0.414861358363668627E-1,(-0.105056666208347189E-1),0.250895935547641537E-2,(-0.481852212924460022E-3),0.542549547083868702E-4,0.578205301871228640E-5;(-0.865350827010216586E0) ...
  ,0.210170633360392454E0,0.113482182369897384E0,(-0.463418900401036875E0),(-0.545468301888350230E0),0.334946010316338282E1,(-0.321649815057877851E1),(-0.141694640067693926E2);(-0.869243889786349909E0), ...
  0.205449185353609747E0,0.138144361797270476E0,(-0.410458595800943704E0),(-0.872238524994977748E0),0.358279073626038449E1,(-0.673975291389731011E0),(-0.250210667636039334E2);(-0.872987559806219154E0), ...
  0.199952311043101368E0,0.158869932914119009E0,(-0.333853289583463593E0),(-0.120032842765269856E1),0.345101508848130838E1,0.334602731899458659E1,(-0.372677817453967673E2);(-0.876703275845851250E0), ...
  0.193632022699006540E0,0.174985898802876102E0,(-0.231549010267667942E0),(-0.149874770413960832E1),0.276487905367067645E1,0.903905771565000390E1,(-0.482179645447263216E2);(-0.880469375815744238E0), ...
  0.186474069505409410E0,0.185027767842526921E0,(-0.103592025162870544E0),(-0.170587884184476203E1),0.127657681346917712E1,0.160744910601467019E2,(-0.509714159048002472E2);(-0.883854825887894667E0), ...
  0.179566534463663556E0,0.187224021691210764E0,0.247518784959613357E-1,(-0.173056964784447126E1),(-0.852153696454206963E0),0.219259682687458339E2,(-0.355607955906770344E2);(-0.886688996862613583E0), ...
  0.173603666107069592E0,0.183414880505823585E0,0.131779221397406580E0,(-0.157321179045427312E1),(-0.310594791108960746E1),0.241456468641801769E2,(-0.719037364994877997E-1);(-0.889271851445130530E0), ...
  0.168169734309373702E0,0.175408094548609193E0,0.218152364768845714E0,(-0.125699439999989193E1),(-0.522159932080078734E1),0.214482675416714008E2,0.540853446041103629E2;(-0.891719539026753781E0), ...
  0.163143149560207167E0,0.164272276066392868E0,0.279775972203414387E0,(-0.808597137643971600E0),(-0.678403781518798698E1),0.125704840529395902E2,0.117854987570814842E3;(-0.894193884406349969E0), ...
  0.158292776892582581E0,0.150453069401885935E0,0.312829842037723342E0,(-0.258378525617612346E0),(-0.726385165716724418E1),(-0.321531047690838051E1),0.169171886834271973E3;(-0.896497348917896437E0), ...
  0.154059562824942281E0,0.136499129504299162E0,0.312760014633940878E0,0.247496051389389952E0,(-0.618679309249265325E1),(-0.210357675147097217E2),0.163721231952616372E3;(-0.898463292853346168E0), ...
  0.150690053047567432E0,0.124761687884795023E0,0.290737746154948438E0,0.582714421655643428E0,(-0.404869537638519262E1),(-0.329350603994848659E2),0.880678120450128281E2;(-0.900303576958432291E0), ...
  0.147740146971075524E0,0.114599146211019929E0,0.257119396876450049E0,0.753146601265499835E0,(-0.144767456326063424E1),(-0.353342738317203276E2),(-0.379834621397131108E2);(-0.902213113631502430E0), ...
  0.144872926873178605E0,0.105317527256148745E0,0.216836625879271049E0,0.762795090852556035E0,0.102327691715492009E1,(-0.255715162492313366E2),(-0.166774756399789705E3);(-0.904028238943482171E0), ...
  0.142307978607487995E0,0.977957532349014234E-1,0.180771713096501457E0,0.649686436561131354E0,0.233573621916624856E1,(-0.848969211090645081E1),(-0.198646249595773847E3);(-0.905640334143535231E0), ...
  0.140142129712329739E0,0.920765969925631772E-1,0.154287858865987269E0,0.509511551385642167E0,0.242225274728328758E1,0.479269585819524907E1,(-0.120205641680722374E3);(-0.907352008805240459E0), ...
  0.137941556803719393E0,0.867985044102688668E-1,0.132609714119970420E0,0.377283559278139952E0,0.181294806657186332E1,0.977938392572631120E1,(-0.877009158721999429E0);(-0.909051540426399237E0), ...
  0.135844833029409256E0,0.821767937272574609E-1,0.116302827734031975E0,0.286277249472785615E0,0.115448497453921218E1,0.714523269156378099E1,0.458670995989005102E2;(-0.910853370402378287E0), ...
  0.133707220746464307E0,0.777943515283471900E-1,0.102764313308201307E0,0.224633918604363005E0,0.737698904935371316E0,0.361741828350979498E1,0.259634690308012503E2;(-0.912857774566247332E0), ...
  0.131422057068545261E0,0.734176476348524992E-1,0.906541955975024603E-1,0.178739527436893490E0,0.501345989079749518E0,0.190771281645735793E1,0.961949740330982433E1;(-0.915032860150628949E0), ...
  0.129042174106018939E0,0.691528536570996260E-1,0.799559471167824103E-1,0.143541070164189549E0,0.354255601277510939E0,0.113445485789441065E1,0.452667523413113249E1;(-0.917581957708222945E0), ...
  0.126372490933435377E0,0.646818860866017947E-1,0.697475798243648443E-1,0.113888108868254716E0,0.248773112373962970E0,0.684054425424293083E0,0.226629493443993194E1];
[0,(1/2),0,(-0.624999999999939512E-1),(-0.752768189698386004E-11),0.104166681137035122E-1,(-0.917262497344728519E-7),(-0.178836983862703116E-2);0.398552556139271908E-1,0.498808879421125936E0,( ...
  -0.149048456166741016E-1),(-0.618396273863123733E-1),0.412320353123484467E-2,0.101789706250281046E-1,(-0.985969407446428193E-3),(-0.171930078459782496E-2);0.782586807987569988E-1, ...
  0.495409821890366050E0,(-0.290377866328872944E-1),(-0.599701104754791424E-1),0.793821761505617043E-2,0.951323986537987802E-2,(-0.186826991263433174E-2),(-0.152577035235930056E-2); ...
  0.114636564530285166E0,0.490158354609336992E0,(-0.420179453345137170E-1),(-0.571253006144370995E-1),0.112751898122455560E-1,0.852087215655875829E-2,(-0.258731071900859751E-2),( ...
  -0.124464192741798076E-2);0.150162627684215385E0,0.483131367060118407E0,(-0.541326644261622631E-1),(-0.534013186982596828E-1),0.141618119981282672E-1,0.726038536106606571E-2,(-0.313751279590101884E-2) ...
  ,(-0.901276102822231043E-3);0.185292241036174213E0,0.474350483296382098E0,(-0.654003922826430737E-1),(-0.488806986580382458E-1),0.165602868342878171E-1,0.579091464894003947E-2,( ...
  -0.350367816410136167E-2),(-0.522103086736221094E-3);0.220242299563301916E0,0.463822240257790607E0,(-0.757496427490024222E-1),(-0.436548034722921031E-1),0.184191302362482152E-1, ...
  0.417858781397161924E-2,(-0.367439072098832732E-2),(-0.135313247709482980E-3);0.255200028184896430E0,0.451523296157515121E0,(-0.850891516907348112E-1),(-0.378178922659578839E-1), ...
  0.196923804976384019E-1,0.249314430674052341E-2,(-0.364671172003518874E-2),0.231217287599283121E-3;0.290366713088693461E0,0.437391091506473900E0,(-0.933141807982440522E-1),(-0.314666502012329186E-1), ...
  0.203415144905325806E-1,0.806785744358245141E-3,(-0.342669628295227326E-2),0.551308071628655899E-3;0.326095959731791689E0,0.421256461285941100E0,(-0.100321416373329245E0),(-0.246793436444130406E-1), ...
  0.203336986575759146E-1,(-0.811073501284885800E-3),(-0.302788837664634266E-2),0.802275132326589791E-3;0.363227066790239496E0,0.402632918710626546E0,(-0.106011564449966823E0),(-0.174579168053710469E-1) ...
  ,0.196210850957907777E-1,(-0.230103971815135627E-2),(-0.246439108019035282E-2),0.965254348938022603E-3;0.402280056233655883E0,0.381061881765104176E0,(-0.110087680877143587E0),( ...
  -0.991108165182684940E-2),0.181417354874795349E-1,(-0.356744376047101275E-2),(-0.176499998956435127E-2),0.101983810513238083E-2;0.440097026628012083E0,0.358305462233123224E0,(-0.112030975308757080E0), ...
  (-0.289560259515710372E-2),0.160782185264038634E-1,(-0.442875408751359607E-2),(-0.104943535561405076E-2),0.962819604667460599E-3;0.473927631124295788E0,0.336459643743233007E0,(-0.112004164820525944E0) ...
  ,0.293120837306799215E-2,0.138027269803572522E-1,(-0.485735854652476981E-2),(-0.433835472245302864E-3),0.833891584813938098E-3;0.504697973581076822E0,0.315433924795817834E0,(-0.110477060552718692E0), ...
  0.770566530139707282E-2,0.114753455675205336E-1,(-0.495630745035795880E-2),0.644689887858114305E-4,0.670159031120336427E-3;0.533416225890437905E0,0.294878114303888880E0,(-0.107732388441335074E0), ...
  0.115898049357723300E-1,0.916981104426600022E-2,(-0.480548943675025928E-2),0.448296895021392135E-3,0.494819211486441776E-3;0.560517339379569760E0,0.274714931532533536E0,(-0.103965432545157084E0), ...
  0.146546516854609646E-1,0.695756418385712234E-2,(-0.446576799774834555E-2),0.720582359582920213E-3,0.325172317507561537E-3;0.586275592291129880E0,0.254924183998082802E0,(-0.993329007290516098E-1), ...
  0.169538585954551236E-1,0.489673129073727498E-2,(-0.399043010084449034E-2),0.888858602068411829E-3,0.173476263936455372E-3;0.610903142310551645E0,0.235495244813520817E0,(-0.939652603730771728E-1), ...
  0.185375384821430554E-1,0.303177290937889213E-2,(-0.342683728117921163E-2),0.964625943900294431E-3,0.475345868288391011E-4;0.634604409141019883E0,0.216397091691486943E0,(-0.879664020800238713E-1), ...
  0.194552354067539082E-1,0.139391040246200275E-2,(-0.281615040038010886E-2),0.962015695046758315E-3,(-0.487068546345582212E-4);0.657662392615029993E0,0.197514992026432492E0,(-0.813926203885472582E-1), ...
  0.197531453896944375E-1,0.318707507636629821E-6,(-0.219156875963823753E-2),0.896068361581568317E-3,(-0.114494736293437198E-3);0.680658187906161475E0,0.178472917443584095E0,(-0.741754314420669965E-1), ...
  0.194567839668776305E-1,(-0.114772141294299052E-2),(-0.157394918192116556E-2),0.779881882291446517E-3,(-0.151269932088884353E-3);0.703507869660555169E0,0.159445427419521657E0,( ...
  -0.664348523574705378E-1),0.185837029903785690E-1,(-0.201220081503446078E-2),(-0.100069780763556247E-2),0.629988404472494743E-3,(-0.160754956446699366E-3);0.724410424616661513E0, ...
  0.142051891645210707E0,(-0.589516955056587608E-1),0.173051761019209338E-1,(-0.253958994000569871E-2),(-0.540102864015979441E-3),0.478062684530274543E-3,(-0.149406947271185187E-3); ...
  0.742447868313190996E0,0.127143204781421409E0,(-0.522718638526649970E-1),0.158665573374474878E-1,(-0.278508637075559246E-2),(-0.209421705712030129E-3),0.347246395628939975E-3,( ...
  -0.128168094226514581E-3);0.758335591810822155E0,0.114163195251364532E0,(-0.462907807031625342E-1),0.143761905709607907E-1,(-0.284245628041930968E-2),0.212074611434212286E-4,0.240042783665686305E-3,( ...
  -0.104040247564539082E-3);0.772677018724104436E0,0.102630993524464677E0,(-0.408760558135886133E-1),0.128840488733594457E-1,(-0.277326732613464674E-2),0.176439136186739759E-3,0.154657442556411300E-3,( ...
  -0.804833081839005604E-4);0.785794547459177045E0,0.922902289344596763E-1,(-0.359649552512535498E-1),0.114280230522079438E-1,(-0.261874146220806084E-2),0.273440951589160951E-3,0.889855785166999986E-4,( ...
  -0.593466353382292842E-4);0.797897849957815596E0,0.829722716730431242E-1,(-0.315169998622193848E-1),0.100358367933844448E-1,(-0.240926125744492156E-2),0.325836154273659337E-3,0.405057061131770206E-4,( ...
  -0.414977822703107037E-4);0.809873186829898312E0,(-0.665428337904364076E0),(-0.206927526719237772E0),0.943017106807903520E0,0.982881557044423939E0,(-0.683446952622129117E1),0.716388437561450346E1, ...
  0.285416336460810449E2;0.822255922846380394E0,(-0.656716127495576049E0),(-0.257377976188469540E0),0.844620655612726464E0,0.165287321161163390E1,(-0.737997697108866139E1),0.200115968677791146E1, ...
  0.510430281136188281E2;0.834288116346279736E0,(-0.646394268112303667E0),(-0.300314816161729521E0),0.697366777715248748E0,0.233192381905563551E1,(-0.717974866431482481E1),(-0.622157026631474138E1), ...
  0.763531832999662042E2;0.846378424615989984E0,(-0.634378168444483409E0),(-0.334353862785579937E0),0.496903759646163950E0,0.295702136841982495E1,(-0.584541625481141905E1),(-0.178949412099086937E2), ...
  0.989194621347713727E2;0.858812456318866429E0,(-0.620640301563995831E0),(-0.356517453944476488E0),0.242943600295117752E0,0.340292951936777774E1,(-0.286904388999328465E1),(-0.323351648409079922E2), ...
  0.104694786791664460E3;0.870169321203950431E0,(-0.607289313062400168E0),(-0.362892876792852020E0),(-0.141370460886731134E-1),0.347996756698014904E1,0.142849131003367219E1,(-0.443805810151699070E2), ...
  0.735619442992747046E2;0.879823605493063282E0,(-0.595707580907024846E0),(-0.356937657771297958E0),(-0.230061985273006072E0),0.318501059874970632E1,0.599912398363004451E1,(-0.490566877532575996E2), ...
  0.184983560841451473E1;0.888747844603027586E0,(-0.585116402722215422E0),(-0.342383478592893404E0),(-0.405571899430677294E0),0.256585277840584608E1,0.103066016311607752E2,(-0.437910472219700346E2),( ...
  -0.107480284164371966E3);0.897321424452189606E0,(-0.575293035496219863E0),(-0.321402642577365158E0),(-0.532124667075567356E0),0.167617202901369457E1,0.135122933474572178E2,(-0.260466172328503524E2),( ...
  -0.236206065543747721E3);0.906106161513912072E0,(-0.565793607496063409E0),(-0.294946088527079266E0),(-0.601932736238515547E0),0.576909138981285906E0,0.145510835977595779E2,0.564824223569220230E1,( ...
  -0.340121480170520848E3);0.914390728188583540E0,(-0.557488301829143194E0),(-0.268011028129129244E0),(-0.605244755981879352E0),(-0.438891096957268046E0),0.124564988372198056E2,0.415284788721092551E2,( ...
  -0.330246069566318011E3);0.921541306015304346E0,(-0.550868335455350509E0),(-0.245266986457358577E0),(-0.563937383080342887E0),(-0.111581292570163348E1),0.821191058456718706E1,0.656038483566957341E2,( ...
  -0.179360783694490082E3);0.928299883328578131E0,(-0.545065880595448502E0),(-0.225547820172572735E0),(-0.498967072963955403E0),(-0.146413695017628814E1),0.301937192954395872E1,0.706732087973833686E2, ...
  0.730448715608836592E2;0.935377490685908242E0,(-0.539419293310357930E0),(-0.207542262351951721E0),(-0.420394363442462239E0),(-0.149130518771530371E1),(-0.193227042570790264E1),0.513675987233001447E2, ...
  0.331773237630622228E3;0.942164677030180349E0,(-0.534361493214235053E0),(-0.192971678429698587E0),(-0.349826109113731795E0),(-0.127167077996888761E1),(-0.457862778325816375E1),0.172954510078525372E2, ...
  0.397011558991793501E3;0.948240276528200554E0,(-0.530085187316743916E0),(-0.181914297927741452E0),(-0.298011838102994757E0),(-0.996001971120573853E0),(-0.477306905830665753E1),( ...
  -0.929137713155092636E1),0.241142541264423817E3;0.954739359781199746E0,(-0.525734530196342628E0),(-0.171728583443954806E0),(-0.255691272804475145E0),(-0.735330742184929194E0),(-0.357334830057260795E1) ...
  ,(-0.193455846600769893E2),0.275892337966810886E1;0.961240819941656310E0,(-0.521583213923344000E0),(-0.162823207178160379E0),(-0.223959193697969666E0),(-0.556151457895272257E0),( ...
  -0.226928840765578395E1),(-0.141531282772500912E2),(-0.910388439611482389E2);0.968185819304771078E0,(-0.517344425768720361E0),(-0.154388428597966998E0),(-0.197694332660403769E0),( ...
  -0.435199861297095646E0),(-0.144458555122246388E1),(-0.714529260027782796E1),(-0.515646994382004290E2);0.975974181959724463E0,(-0.512805210992120470E0),(-0.145972118823086936E0),( ...
  -0.174260151668745794E0),(-0.345471641982822688E0),(-0.978574717377344602E0),(-0.375423662331726548E1),(-0.190405445258837892E2);0.984499694408579218E0,(-0.508068596186012462E0),( ...
  -0.137776671309925339E0),(-0.153603822735611525E0),(-0.276861605039234512E0),(-0.689537562370305807E0),(-0.222578393719901838E1),(-0.893322597005477198E1);0.994588434186297253E0,( ...
  -0.502743121302737814E0),(-0.129189492334971191E0),(-0.133933846609077161E0),(-0.219221551516641177E0),(-0.482880294246942834E0),(-0.133810464828781443E1),(-0.445955548956127013E1)];
[0,1,0,(-3/8),0.109188995499186608E-11,0.104166665296848485E0,0.274977775251019606E-6,(-0.250844256683920160E-1),0.554229075445027057E-3;0.795838622182803433E-1,0.992861688014737726E0,( ...
  -0.892185130480711933E-1),(-0.368406371564524584E0),0.411057966471632037E-1,0.100845174339058329E0,(-0.137611551988019208E-1),(-0.237869598270102964E-1),0.389538099436418170E-2;0.155557803408837237E0, ...
  0.972584017013087893E0,(-0.172642670783284147E0),(-0.349850316741758623E0),0.784414516864382841E-1,0.916104730197257704E-1,(-0.257878172006987427E-1),(-0.203427336463664238E-1), ...
  0.713750353911837583E-2;0.226253505712736364E0,0.941526357711030550E0,(-0.247177642174251037E0),(-0.321933621397595356E0),0.109867389787830076E0,0.780390243421346073E-1,(-0.350881028422383235E-1),( ...
  -0.154245569025992443E-1),0.936205337205314614E-2;0.293527683597711882E0,0.900485859839519718E0,(-0.313862976288299870E0),(-0.285991730071693437E0),0.135349835098145884E0,0.611605487067908339E-1,( ...
  -0.415041739694837100E-1),(-0.957271702291378008E-2),0.105025547502608817E-1;0.357787696879363034E0,0.850042035890923360E0,(-0.372209104344453432E0),(-0.243320566845519010E0),0.154321833606884738E0, ...
  0.420447805405908610E-1,(-0.448222614766757034E-1),(-0.334728515942090175E-2),0.105348058161835068E-1;0.418943337179410805E0,0.790804090516322232E0,(-0.421291016658063863E0),(-0.195381149050301622E0), ...
  0.166224416206987293E0,0.218630082329473819E-1,(-0.449684179567935609E-1),0.267724344667057998E-2,0.953151365422456244E-2;0.476792558078836469E0,0.723344350393690555E0,(-0.460159658443638603E0),( ...
  -0.143729530404002331E0),0.170702413608686562E0,0.181601026630753046E-2,(-0.420685971649921139E-1),0.796456023924229145E-2,0.766783436126002671E-2;0.531067758123930683E0,0.648183362139535743E0,( ...
  -0.487874444285524869E0),(-0.900075527256479585E-1),0.167630005119359416E0,(-0.169063178451478445E-1),(-0.364472184945092316E-1),0.120544558677966156E-1,0.520445680281668187E-2;0.581577420113030467E0, ...
  0.565509833581126427E0,(-0.503500945204937259E0),(-0.357870233319246912E-1),0.157070842534130763E0,(-0.332014248881285763E-1),(-0.285879485836709129E-1),0.145994374150185938E-1, ...
  0.245519219779905288E-2;0.628434673673907637E0,0.474337392104957719E0,(-0.505791885695646841E0),0.177518797903299818E-1,0.139011216100788670E0,(-0.460953133647545671E-1),(-0.190254383349910181E-1), ...
  0.153613749669203184E-1,(-0.255777872481638212E-3);0.670709337515912024E0,0.374591769158410800E0,(-0.492684388928663250E0),0.682590258941937867E-1,0.113738473332111037E0,(-0.543135955263379268E-1),( ...
  -0.860404511335107522E-2),0.141887198483501712E-1,(-0.253266819637713278E-2);0.703955132797613703E0,0.276401168192606286E0,(-0.465190688016712187E0),0.108980514757046980E0,0.851203367711478576E-1,( ...
  -0.566509418537119615E-1),0.662991352498392773E-3,0.114941037211960366E-1,(-0.389008270512522622E-2);0.726568918962128997E0,0.189183740202836568E0,(-0.429027235686860370E0),0.136812621766128487E0, ...
  0.579756048375132901E-1,(-0.541831052053014453E-1),0.741390976463188647E-2,0.825341645372350228E-2,(-0.430194396065465049E-2);0.740725056346424546E0,0.112007191497134550E0,(-0.387623218475210027E0), ...
  0.154022968082721692E0,0.336089303964689079E-1,(-0.486404576806654754E-1),0.118031080267745482E-1,0.506336384123092241E-2,(-0.406030220543789652E-2);0.747964235753955452E0,0.432255118320229453E-1,( ...
  -0.342736237567443335E0),0.162576451658071595E0,0.124123774364584938E-1,(-0.412165237234423121E-1),0.141767667622579036E-1,0.223085619947347641E-2,(-0.341562019335082177E-2);0.749114410713419402E0,( ...
  -0.175724713233317210E-1),(-0.295977339108583299E0),0.163817640830525235E0,(-0.522754229164062668E-2),(-0.328521182965317502E-1),0.148627061111548886E-1,(-0.541664727537532173E-4),( ...
  -0.257082465303344385E-2);0.744753033279814650E0,(-0.705475596222747253E-1),(-0.248741099455497472E0),0.158945725770693292E0,(-0.191157393244179156E-1),(-0.243085832895527308E-1), ...
  0.142199396194301114E-1,(-0.170780251514202533E-2),(-0.168750843797660108E-2);0.735313114271660971E0,(-0.115805972378285878E0),(-0.202215558802385239E0),0.149090989496010211E0,( ...
  -0.292458044749763059E-1),(-0.161832790560048563E-1),0.126177494705210767E-1,(-0.273043165215595977E-2),(-0.883128205483578172E-3);0.721101100789166898E0,(-0.153467571579917642E0),( ...
  -0.157373276308035131E0),0.135311137524303972E0,(-0.357700769112175945E-1),(-0.891627370647598403E-2),0.104088454903906956E-1,(-0.318352898236060588E-2),(-0.231462333130778729E-3); ...
  0.702218410758784849E0,(-0.183714909689010801E0),(-0.114887400427511340E0),0.118523393868181146E0,(-0.389555015237887906E-1),(-0.279331122270377662E-2),0.790203509130622093E-2,( ...
  -0.316579844886838350E-2),0.233104305663164317E-3;0.678215332265644083E0,(-0.206802610554785635E0),(-0.748884211215600919E-1),0.992948616599767195E-1,(-0.390876037131532375E-1), ...
  0.204648836861571214E-2,0.533131072082626608E-2,(-0.278551893534279116E-2),0.506771098630727374E-3;0.649079694046669138E0,(-0.222002845785103356E0),(-0.387845613766988726E-1),0.787366840262320010E-1,( ...
  -0.364660896035967390E-1),0.537998598899278538E-2,0.297434493150107472E-2,(-0.217099059922196595E-2),0.598964168846171947E-3;0.617711738284921797E0,(-0.228598355120356837E0),(-0.100524406918891439E-1) ...
  ,0.596575021497517931E-1,(-0.320598807556692375E-1),0.706810909872885938E-2,0.118493662249551600E-2,(-0.152086086441060317E-2),0.555975224785368767E-3;0.586996242062410378E0,(-0.228372301142625573E0), ...
  0.106712697399540257E-1,0.437664757790166821E-1,(-0.271149884162058360E-1),0.751946399349105694E-2,0.218733332636230213E-4,(-0.976238480165290876E-3),0.453513923979927286E-3;0.557168992278572588E0,( ...
  -0.223514285113289393E0),0.253241781653998145E-1,0.307671688435441939E-1,(-0.222221353984083886E-1),0.723535477203611211E-2,(-0.676032396026108252E-3),(-0.558481322910508936E-3), ...
  0.338920933508372850E-3;0.528073603873753367E0,(-0.215382213393799458E0),0.353750593294625865E-1,0.202264478028618554E-1,(-0.176469198284080812E-1),0.653320750960903438E-2,(-0.104445125100185075E-2),( ...
  -0.256602479229852344E-3),0.233674142929471109E-3;0.499738031877073784E0,(-0.204909440716034253E0),0.417833770108666971E-1,0.118477475158734483E-1,(-0.135467170120741748E-1),0.562006195362933196E-2,( ...
  -0.118263367307237772E-2),(-0.534168700249016554E-4),0.147076969313619367E-3;0.472229142962228924E0,(-0.192807736994663979E0),0.452861594387038143E-1,0.536666474924493004E-2,(-0.100017547867551904E-1) ...
  ,0.463675474800711591E-2,(-0.116775951792592117E-2),0.709540925351595801E-4,0.815750971545304932E-4;0.443843117946304104E0,0.160689949937273085E1,(-0.105927739781994926E1),(-0.828044134105229357E1), ...
  0.149300306656238915E2,0.395792270060927839E2,(-0.219325941746507605E3),0.308935333346701844E3,0.775424646829201551E3;0.413435625037939142E0,0.163749685456656433E1,(-0.565680846268626235E0),( ...
  -0.922998531756433930E1),0.100073441696587877E2,0.661340269740649185E2,(-0.249187130787104879E3),0.117198130149579228E3,0.185441292885185261E4;0.383068569700178872E0,0.164873589761692690E1,( ...
  -0.385689979302448005E-1),(-0.971201040783249886E1),0.261907105142776872E1,0.938167342206069814E2,(-0.242504582222292094E3),(-0.259552621157738170E3),0.329738182936235920E4;0.351997798420947386E0, ...
  0.163980259743251168E1,0.510263580918957714E0,(-0.954447674834854435E1),(-0.743889426731019773E1),0.117959549096372760E3,(-0.170171399714595993E3),(-0.878915800514512727E3),0.487131445543910635E4; ...
  0.319782083785235379E0,0.160866858630151000E1,0.105054307028297519E1,(-0.847103903146845765E1),(-0.198320847358871542E2),0.128695467224805395E3,0.950306086955069211E1,(-0.173593345672640961E4), ...
  0.566867722594280087E4;0.290436998922798002E0,0.156168672043337410E1,0.147185469569761074E1,(-0.657238038138943572E1),(-0.312556157420604357E2),0.113217695415394047E3,0.285858099421118527E3,( ...
  -0.248525706074421113E4),0.385118718078216186E4;0.265774833598219198E0,0.150991050969925088E1,0.173567697519173818E1,(-0.430360162154383687E1),(-0.388626263761232982E2),0.715645569386638117E2, ...
  0.583010679154054354E3,(-0.267747961916275924E4),(-0.159064647473594524E4);0.243360090896239624E0,0.145503929555342399E1,0.187558714003448467E1,(-0.183530637582102806E1),(-0.419603162203721748E2), ...
  0.650042049036094370E1,0.839295600389179193E3,(-0.196064779766375181E4),(-0.109197135229397180E5);0.222270862748668685E0,0.139894104230968684E1,0.190238742100632578E1,0.602311112256330843E0,( ...
  -0.395125426887683618E2),(-0.744741281274543355E2),0.953227119338942256E3,0.799452037024788708E1,(-0.224395266190980427E5);0.201175365114969681E0,0.134133114376490038E1,0.182185967893945929E1, ...
  0.279108219074346721E1,(-0.304845048488466700E2),(-0.157424140313653084E3),0.781769053175334700E3,0.335149341814695166E4,(-0.302111032488098132E5);0.181773736887434427E0,0.128975109111211702E1, ...
  0.166407950070772377E1,0.420388759472027682E1,(-0.167966334730973539E2),(-0.205945609701263455E3),0.256976983119236891E3,0.663748717691353224E4,(-0.218064931925957024E5);0.165398247140012000E0, ...
  0.124901952245779674E1,0.148903463400297485E1,0.472372053463785969E1,(-0.340162683563709025E1),(-0.200597575123749353E3),(-0.410564090281561793E3),0.761723000445771227E4,0.608141642515850155E4; ...
  0.150209476774239508E0,0.121444488096644619E1,0.131470792646821405E1,0.460775777031394117E1,0.755273640873659316E1,(-0.147400419387792765E3),(-0.991497416945030474E3),0.514491241127151766E4, ...
  0.442895594836439305E5;0.134570857787745791E0,0.118239120556347730E1,0.114480089907719994E1,0.401053627295187582E1,0.143407947411882763E2,(-0.578264825685637357E2),(-0.120068552862366930E4),( ...
  -0.100835640911477205E4),0.670229334423541113E5;0.119797864402112210E0,0.115524716625898561E1,0.100716800528940255E1,0.323933324890963205E1,0.153057801737418038E2,0.225808974246657434E2,( ...
  -0.832486510203287957E3),(-0.675841712505621621E4),0.367547430289373381E5;0.106736303596193973E0,0.113342976840223916E1,0.907677400512462579E0,0.259027587652668456E1,0.127734522857468679E2, ...
  0.591480591714204681E2,(-0.227214535322995330E3),(-0.734700307208191262E4),(-0.249145730799660842E5);0.929151510252532323E-1,0.111216960925675176E1,0.822482335821888121E0,0.205397603193172106E1, ...
  0.904023512156205812E1,0.564271993444733889E2,0.232481686857213415E3,(-0.275975772835425806E4),(-0.558439056123235058E5);0.792296370476952895E-1,0.109263283421606366E1,0.753352510654611643E0, ...
  0.168167765449804865E1,0.617274738198265901E1,0.355940522627125963E2,0.267564714469928657E3,0.123783108340830458E4,(-0.177824774713842760E5);0.647519274966606369E-1,0.107333601901648584E1, ...
  0.691785242376913911E0,0.140404908104379373E1,0.438565661486083484E1,0.197730892320854703E2,0.130073705207384072E3,0.116990831263984805E4,0.971304695699174309E4;0.486750627994973755E-1, ...
  0.105332153684070834E1,0.633510380281155984E0,0.117672868398071716E1,0.322820727818728126E1,0.119153812805148714E2,0.574347163657633563E2,0.365999110032306988E3,0.310244579923098961E4; ...
  0.312527966714801753E-1,0.103308729055420247E1,0.579452556665423080E0,0.989745325783654818E0,0.242697407680094100E1,0.771411022443057463E1,0.304488476119510583E2,0.146907248920638255E3, ...
  0.844269608526141512E3;0.108529026926124341E-1,0.101106396581930706E1,0.525431808550799004E0,0.822532774717782700E0,0.180589551228619770E1,0.500157297109078878E1,0.166794698479250369E2, ...
  0.655691890391038597E2,0.295166054498520927E3]
};

polorder = size(funcs{1},2);
disp(['const int numPolynomes = ' num2str(numel(origins))        ';']);
disp(['const int polyorder = ' num2str(polorder) ';']);
disp(['const int numTableElements= numPolynomes*polyorder;']);
disp(['const int cutoffbinidx = ' num2str(cutoffbinidx) ';']);

tmp1 = reshape(typecast( single(limits(:)) ,'uint32'),1,[]);
str1 = sprintf('0x%08X, ', tmp1);
fprintf('ALIGN_VEC const long tablelimits[(numPolynomes+1)] = { %s };\n', str1(1:end-2))
tmp2 = reshape(typecast( origins(:)     ,'uint32'),2,[]);
str2 = sprintf('0x%08X,0x%08X  , ', [tmp2]);
fprintf('ALIGN_VEC const long origins[(numPolynomes)*2] = { %s };\n', str2(1:end-2))

for k=1:numel(funcs)
    tmp1 = reshape(typecast( reshape(funcs{k}(:,end:-1:end-polorder+1)',[],1) ,'uint32'),2,[]);
    str1 = sprintf('0x%08X,0x%08X  , ', [tmp1]);
    fprintf('ALIGN_VEC const long fun%d[(numTableElements)*2] = { %s };\n', k, str1(1:end-2))
    if size(funcs{k},2)>polorder
        tmp1 = reshape(typecast( reshape(funcs{k}(:,end-polorder:-1:1)',[],1) ,'uint32'),2,[]);
        str1 = sprintf('0x%08X,0x%08X  , ', [tmp1]);
        fprintf('const int fun%d_numextra = %d;\n',k, size(funcs{k},2)-polorder);
        fprintf('ALIGN_VEC const long fun%d_extra[(fun%d_numextra * numPolynomes)*2] = { %s };\n', k, k, str1(1:end-2))
    end;
end;

function testmex()
%%
% bx = linspace(-100,100,1000000)';
bx = randn(1000)*270;
tic;
[a]=rice_besseli_funs_c(bx);
toc,tic;
tic;
[a,b]=rice_besseli_funs_c(bx);
toc,tic;
tic;
[a,b,c]=rice_besseli_funs_c(bx);
toc,tic;
besX = besseli(0,bx,1);
toc,tic
at = log(besX);
toc,tic
bt =  besseli(1,bx,1)./besX;
toc,tic;
ct =  (besseli(2,bx,1)./besX + 1 - 2 * bt.^2).*bx;
toc
%%
[dummy, reord]= sort(bx(:));
plot(bx(reord),[a(reord)-at(reord) b(reord)-bt(reord) c(reord)-ct(reord)])
%%
sum( abs(reshape(a-at,[],1))./eps(at(:)) >100)/numel(at)
sum( abs(reshape(b-bt,[],1))./eps(bt(:)) >100)/numel(bt)
sum( abs(reshape(c-ct,[],1))./eps(ct(:)) >100)/numel(ct)

[sum(~isfinite(a(:))), sum(~isfinite(at(:)))]
[sum(~isfinite(b(:))), sum(~isfinite(bt(:)))]
[sum(~isfinite(c(:))), sum(~isfinite(ct(:)))]

% NOTE: there are elements with larger difference, but this (mainly) is due to besseli beiing inaccurate
%%
bx = [-1 0 1 1e40 inf nan];
besx=log(besseli(0,bx,1))
[a]=rice_besseli_funs_c(bx)
%% test final rice distribution:
A=linspace(-10,10,1000)';
x = 4;sigma = 2.345;
p0 = ricepdf(x,A,sigma);
[a b c]=logricepdf(x,A,sigma,[false true false]);
plot(A,[a b c log(p0)],(A(1:end-1)+A(2:end))/2,diff([a b])/median(diff(A)))

%%
if 0 
    % some test to validate that we found the correct optimum =>
    % derivatives are correct.
tstbase = truefulltht(:,show_param_index)

rng1=linspace(0,4,100)';
rng2=linspace(-3,4,100)';
rng3=linspace(0,4,100)';
rng4=linspace(.5,2,100)';
L = zeros(numel(rng1),numel(rng2),numel(rng3),numel(rng4));
%%
kstart = 30742460;nblock = 1000;
for kst=kstart : nblock: numel(L)
    ked = min(kst+nblock-1, numel(L));
    k= (kst:ked)';
    [i1,i2,i3,i4] = ind2sub(size(L), k);
    x = [rng1(i1), rng2(i2), rng3(i3), rng4(i4)]';
    L(k) = sum( logricepdf(simdata(:,show_param_index,28)*ones(1,size(x,2)), fun(x(1:3,:)), x(4,:)) ,1);
end;
%%
imagebrowse(L);
end;